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Abstract. In this paper we explore the evolution 
of a PWN while the pulsar is spinning down. An 
MHD approach is used to simulate the evolution of 
a composite remnant. Particular attention is given 
to the adiabatic loss rate and evolution of the neb- 
ular field strength with time. By normalising a two 
component particle injection spectrum (which can 
reproduce the radio and X-ray components) at the 
pulsar wind termination shock to the time dependent 
spindown power, and keeping track with losses since 
pulsar/PWN/SNR birth, we show that the average 
field strength decreases with time as <^^, so that the 
synchrotron flux decreases, whereas the IC gamma- 
ray flux increases, untU most of the spindown power 
has been dumped into the PWN. Eventually adiabatic 
and IC losses wiU also terminate the TeV visibility 
and then eventuaUy the GeV visibUity. 

I. Introduction 

Aharonian et al. fT! discussed eight unidentified VHE 
gamma-ray sources discovered with the H.E.S.S. tele- 
scopes. All are extended objects with angular sizes 
ranging from approximately 3 to 18 arc minutes, lying 
close to the Galactic plane (suggesting they are located 
within the Galaxy). In each case, the spectrum of the 
sources in the TeV energy range can be characterized 
as a power-law with a differential spectral index in 
the range 2.1 to 2.5. The general characteristics of 
these sources (spectra, size, and position) are similar to 
previously identified galactic VHE sources (e.g. pulsar 
wind nebulae PWNe), however since these sources have 
so far no clear counterpart in lower-energy wavebands, 
further multi-wavelength study is required to understand 
the emission mechanisms powering them, and therefore 
follow-up observations with higher-sensitivity X-ray and 
GeV 7-ray telescopes will be beneficial (as stated in 1 1 ].) 

One possibility is that we are dealing with relatively 
old PWN born from Type II supernovae, but still rel- 
atively close to the molecular clouds from which the 
massive progenitor stars were born. This will then also 
explain their proximity to the galactic plane. 

A natural explanation would be that these sources 
were once bright in synchrotron emission, but that the 
field strength decreased with time as the PWN expanded 



with time |2|: The pulsar eventually deposited all its 
spindown power into the nebula and whereas the syn- 
chrotron brightness decreased with time because of field 
decay, the inverse Compton 7-ray flux increases until 
reaching a convergent value, after which it will also 
decay because of continuous adiabatic losses and inverse 
Compton cooling. The 7-ray lifetime of a PWN can 
then be much longer than the apparent radio and X-ray 
lifetimes. 

In this paper we will discuss the results of MHD 
and radiative modelling of evolving PWNe and show 
predicted evolutionary results for the composite SNR 
G21.5-0.9. 

II. The MHD model for composite supernova 

REMNANTS 

Supernova remnant evolution in either uniform or non- 
uniform media have been modelled extensively by e.g. 
0, Q. For eiflier composite SNRs or PWNe in the ISM 
simulations were also presented by e.g.f6l, fl], fE\, 19], 
1 10 1, 1 11]. In this work we use a similar model as used 
in most of the studies above by solving the well known 
Euler equations 

^p + V-(pv) = 0, (1) 

|^(pv) + V-(pvv + PI) = 0, (2) 

ot 



dr2 



+V-^v2v+^, 

7 — 1 2 7 — 1 



= 



(3) 



which describe inviscid flow. Here p is the density, v 
the velocity and P the gas pressure. These equations 
describe the balance of mass, momentum and energy. 
Currently we only consider a one fluid scenario with an 
adiabatic index of 5/3. Although a relativistic description 
is necessary to model PWN evolution correctly, the 
speed of the relativistic material downstream of the 
pulsar wind termination shock is sufficiently smaller 
than c to use a non-relativistic treatment (see also e.g. 
1 6 1). The numerical scheme is discussed in |17| and 
compute solutions to hyperbolic differential equations 
using a wave propagation approach. See also 111 811 for 
more details. The model solves in spherical coordinates 
r and (p, with r ranging from 0.01 pc to 25 pc (2000 
gridpoints) and <j> from 0° to 180° (150 gridpoints). 
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For the initial and boundary conditions of the SNR 
(see also fl^, f6l, lH), JTOl) we assume a spherical 
region, radius r^j, and a high constant density p^j with 
a radially increasing velocity profile 



7 = ^ej-r/ 



(4) 



In this case we take r^j = 0.1 pc while for the density 
we have 

Pej 



Airr 



(5) 



with Mej the ejecta mass. For the velocity we have 



'10^ 

3 Me 



(6) 



V X (v X B) = 



(7) 



To compute the PWN magnetic field we solve 

dB 

'dt 

using a similar scheme as for the fluid part. Note that 
this is not a full MHD solution because the field is 
calculated kinematically from the flow ( |fT9ll . ||201 ) and 
no backreaction on the fluid is considered. More detailed 
MHD calculations were done by e.g 111 & ifTTl . 

III. Constraints from pulsar evolution 

For the PWN we assume that the spin-down luminos- 
ity of the pulsar is given by (assuming a pulsar braking 
index of 3) 



Lit) 



(8) 



where Lq is the initial spindown power and r the 
spindown timescale, which, for a birth period Po and 
present period P, is defined as 



(9) 



IV. The evolution of the plerionic magnetic 

FIELD STRENGTH 

To calculate the multiwavelength (MWL) spectrum 
we need to know the behaviour of the average PWN field 
strength B{t) with time. This quantity was calculated by 
taking the volume averaged field strength between the 
pulsar wind termination shock radius and PWN outer 
radius. 

The calculation of the average field strength starts 
progressively later (in time) with decreasing io- This 
is because of the difficulty in resolving the position 
of the PWN termination shock as Lq decreases. This 
difficulty should be resolvable if we reduce the grid size 
of the calculation, but at the expense of CPU calculation 
time. For example, the PWN termination shock radius 
of G21.5-0.9 is ^ 0.5 arcsec, corresponding to a shock 
radius of 0.01 pc, which is already consistent with the 
minimum assumed grid size. 

Figure [T] shows the behaviour of the average B{t) 
for io ranging between 10'^® to 10^^ erg/s and ISM 



L„=10" erg/s 
L.-IO- erg/s 
L.-IO- erg/s 
L„=10'= erg/s 





Fig. 1 . The average magnetic field strength of the PWN for t = 3000 
y and Lo in the range as indicated. Two values of the ISM density 
was assumed. The slope oc t~^ '^ indicates the approximate pre-reverse 
shock field decay evolution with time. 



densities of 10^^^ glcw? and 10^^^ g/cm"^. A more 
detailed discussion of this will be given elsewhere. 

Prior to the passage of of the reverse shock, we find 
that the field strength decreases as t^^-^, independent of 
the chosen parameters. This is modified by the reverse 
shock, but after passage, the time evolution is expected 
to revert back to this t^^-^ behaviour 

V. ADIABATIC LOSSES 

In this section the evolution of a PWN inside a SNR is 



8Mp 







studied. Model solutions corresponding to Me 
in Equation |5] and spin-down time r = 3000 y and r 
= 300 y in Equation |8] are shown . Different scenarios 
ranging from initial pulsar wind luminosity Lq ~ 10"'^ 
erg/s to = 10"^^ erg/s in Equation [8] are shown. 

The rate of change of the energy of a particle con- 
vected by a pulsar wind expanding at a velocity V is 
given by 



(10) 



Below we will see that this quantity is expected to be 
negative, giving rise to adiabatic losses, unless the PWN 
is sufficiently crushed by the reverse shock, such that 
the term V • < 0, in which case the particles will 
experience adiabatic heating. For practical purposes we 
calculate the average adiabatic energy loss rate over the 
PWN between the termination shock and PWN radii by 
averaging the quantity V • V over volume. The radius 
of the PWN was determined by establishing the position 
where the PWN field strength drops to zero. We scale 
the abovementioned rate of energy change (averaged 
over volume) by multiplying the relative energy loss 
rate {dE/dt)/E with the age t of the PWN to give 
the dimensionless quantity t{dE/dt)/E. The results are 
shown in Figures |2a| and |2b] for spindown timescales of 
T = 300 y and r = 3000 y respectively and PSR/SNR 
parameters discussed above. 

Initially we find that the quantity t{dE/dt)/E is 
negative as a result of expansion, so that the particles 
loose energy due to this process. However, when the 
reverse shock compresses the PWN, we find that the 
quantity V • V becomes negative, in which case the 
particles will start to gain some of their lost energy. 
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7^300 y, Lq^IO^^ erg/s 
t=300 y. lo=1o''° erg/s 
!=300y, Lq^IO'*'' erg/s 




A 




TIME t (y) 

(a) T = 300 yr 



TIME t (y) 

(b) T = 3000 yr 



Fig. 2. The scaled relative energy loss rate t{dE/dt)/E (a dimensionless quantity) due to adiabatic expansion as a function of time t since 
birth. The spindown timescale in this case is t = 300 y (left) and 3000 y (right), whereas the SNR ejecta mass for both cases was 8Mq. The 
spindown power Lq at birth is indicated in the Legend. 



With decreasing Lq we find that this heating process 
starts at progressively later times, since the time when 
the reverse shock encounters the PWN increases with 
such decreasing Lq for the same r. 

It is interesting to note that the quantity \t{dE/dt) /E\ 
is always less than unity except when the reverse shock 
compresses the flow. 

Since the relative adiabatic energy loss rate is nearly 
constant at a value around -0.5 (excluding the time of 
reverse shock passage), the total adiabatic losses of a 
particle injected during birth with initial energy and 
which can survive significant radiation losses would be 



E^Eo[^ 
to 



-0.5 



(11) 



Note that = if to = 0, which implies an inconsistent 
solution, unless to > 0. We find that to < 100 yr but 
we are currently reducing the grid size of the simulation 
and will report on the solution for a convergent value of 
to in a followup paper 

VI. Time dependent evolution of the lepton 

SPECTRUM 

We define N{E,t) as the time dependent differential 
particle spectrum for leptons of energy E — ^rrieC^, 
whereas Tgyn and Tad are the timescales corresponding 
to synchrotron and adiabatic losses respectively. The 
magnetic field strength B{t) (used in Tgyn) is time 
dependent. We then integrate the transport equation 



dN 
~dt 



N 



'syn 



N 

Tad 



Q{t) 



(12) 



between time i = when P = Po, i.e. the pulsar birth 
period and the current epoch at Tsnr assuming a pulsar 
braking index n — 3. From N we calculate the spectral 
energy distributions (SED) in synchrotron and inverse 
Compton as discussed below. 



We adopt the injection spectrum of 1 16| for electrons 
at the pulsar wind shock 

n(F i\=( Qo{t){E/Et)-P^ for E<E, 

' ^ V Qoit)iE/Eb)-P'^ for Eb<E< S„,ax 

(13) 

with Et the intrinsic break energy between the radio 
and X-ray components. A value of pi ^ 1.0 reproduces 
the typical flat radio spectra, whereas p2 ^ 2 would 
reproduce the uncooled spectral indices seen in X-rays 
at the pulsar wind termination shock. 

Following |16J, the energy equation for Q{t) can be 
written in terms of the time dependent spindown power 
L{t) giving 

J Qij,t)EdE = 7^L{t). (14) 

We will assume the conversion efficiency r] of spindown 
power to particles as a free parameter The total injected 
lepton energy over time t since birth is then (assuming 
a constant 77) 

We{t) = / TjL{t)dt = f^AE,ot, (15) 
Jq 

where AE.^t = - f^^)/2 (with O = 2tt/P) is the 
net kinetic rotational energy deposited between birth and 
time t. 

VII. Evolution towards an unidentified 

GAMMA-RAY SOURCE 

In our evolutionary model we will use the young 
composite SNR G2 1.5-0. 9 as an example and follow 
the time evolution of the leptonic spectrum and hence 
MWL intensity. The central pulsar PSR J1833-1034 has 
a period of 61.8 ms and for an expansion age near 1 
kyr II2TI . the spindown timescale r should vary between 
3000 and 3800 y given an inferred birth period Pq 
between ~ 50 and 55 ms. The corresponding initial 
spindown power ranges between Lq = 5 x 10'^^ and 
10^^ erg/s. 

We will use pi = 1 as observed in radio ifTSl while 
for X-rays we would expect that a value of p2 =2 
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corresponding to the pulsar wind termination shock [[14) 
would reproduce the MWL spectrum best. However, an 
average value of p2 = 2.6 seems to fit the data better. 

To reproduce the ratio of energy fluxes between X- 
rays and TeV, we normalise the average field strength to 
25/iG at the present age near 1 kyr 

ISO observations |[l21 of the PWN show that the 
radio spectrum should break around 10^^ Hz. This break 
is either intrinsic or due to radiation losses. We find 
however that this break cannot be due to radiation 
losses since this would imply a too large Crab-like field 
strength, which cannot be reconciled with the observed 
ratio of TeV to X-ray flux. An intrinsic break at energy 
Eb to p2 — 2.6 best reproduces the post spectal break 
data. 

For a birth period of Pq = 50 ms we still need a 
relatively large conversion efficiency of = 0.7 in eqn 
[T4lto reproduce the observed synchrotron and IC spectra 
at the present epoch. The required break energy in eqn 
[T3] is El, = 40 GeV, which we will keep fixed with 
time since we have no theory on the evolution of Ei,. 
The assumed radiation fields for the IC calculation were 
the CMBR, a 25K galactic dust component and starlight 
component corresponding to 1 eV/cm'^. The latter two 
radiation energy densities agree approximately with the 
values found for the inner galaxy region at the location 
of G2 1.5-0.9 by E). 

Assuming a constant 77 with time, but the spindown 
power decreasing as that of a magnetic dipole, and hence 
decreasing particle input with time, we were able to 
calculate the time evolution of (5(7, t) and hence the 
MWL spectrum from which the time evolution of the 
radio. X-ray and TeV fluxes were calculated. The latter 
two are shown in Figure [3] It is clear that the X-ray flux 
decreases with time given the decreasing magnetic field 
strength with time, whereas both the inverse Compton 
GeV and TeV fluxes increases with time, reaching a 
limiting value. The predicted radio. X-ray and TeV 
fluxes agree with the observed fluxes at the present 
epoch. 



VIII. Conclusions 

In this paper we have given the basic ingredients 
which gives the time evolution of the MWL spectrum of 
a PWN. The basic result is the following: Whereas the 
X-ray flux is large during early epochs, the GeV and TeV 
fluxes start at relatively low values. As time progresses 
towards the Vela and post- Vela epoch, the synchrotron 
flux starts to decrease significantly, whereas the IC flux 
uncreases, until reaching a steady state value. Given the 
page limit of this paper, we could not explore the details 
of IC and adiabatic losses which would affect the time 
evolution at epocs ^ 10 kyr This will be discussed in 
a followup paper 

The basic conclusion however remains, as a PWN 
grows older, it can remain bright in IC, whereas the 
GeV/TeV flux remains high. This can continue until IC 
and adiabatic losses, or, breakup and diffusion into the 
ISM finally terminates the gamma-ray lifetime. 
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